module Interpolation_module
    implicit none
contains
    !*******************************************************************************
    ! Lagrange Interpolation
    ! xData - x-coordinates
    ! yData - y-coordinates
    ! x0    - x coordiante at which f(x) is to be found out
    !*******************************************************************************

    real function phyLagrangeInterpolation(xdata, ydata, x0)
        real, allocatable, intent(in)	::	xdata(:), ydata(:)
        real, intent(in)				::	x0
        real, allocatable               ::	w(:)
        real		              		:: 	x, y
        integer				            ::	i, j, N
        N = size(xdata, dim=1)
        allocate(w(N))
        do j = 1,N
             w(j) = 1
            do i = 1,N
                if(i == j) cycle
                w(j) = w(j)*(x0 - xdata(i))/(xdata(j) - xdata(i))
            end do
            w(j) = ydata(j)*w(j)
        end do
        phyLagrangeInterpolation = sum(w)
    end function phyLagrangeInterpolation

    !*******************************************************************************
    ! Newton's Divided Difference Method(Not working)
    ! xData - x-coordinates
    ! yData - y-coordinates
    ! x0    - x coordiante at which f(x) is to be found out
    !*******************************************************************************

    real function phyNewtonInterpolation(xdata, ydata, x0)
        IMPLICIT NONE
        real, allocatable, intent(in)   ::  xdata(:), ydata(:)
        real, intent(in)                ::  x0
        real, allocatable               ::  D(:,:),xdifference(:)
        integer                         ::  i, j, N
        N = size(xdata, dim=1)
        allocate(D(N,N))
        allocate(xdifference(N))
        !initializing the first column of coeff matrix with the y data itself
        D(:,1) = ydata(:)
        !generating the rest of the coeff matrix
        do j = 2, N
            do i = 1, N-j+1 
                D(i,j) = (D(i+1,j-1) - D(i,j-1))/(xdata(i+j-1)-xdata(i))
            end do
        end do
        do i = 1,N
            print*,D(i,:)
        end do
        print*,""
        xdifference(1) = x0 - xdata(1)
        do i = 2, N
                xdifference(i) = xdifference(i-1)*(x0 - xdata(i))
                print*,xdifference(i)
        end do
        phyNewtonInterpolation = D(1,1) + DOT_PRODUCT(D(1,2:), xdifference(:))
    end function phyNewtonInterpolation
end module Interpolation_module
